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Abstract 

The  well-known  Stormer-Cowell  class  of  linear  k-step  methods 
for  the  solution  of  second  order  initial  value  problems  suffer 
from  orbital  instability.  The  solution  of  a  problem  describing  a 
uniform  motion  in  a  circular  orbit,  spirals  inward.  Several 
modifications  were  suggested,  some  require  the  a-priori  knowledge 
of  the  frequency.  Here  we  develop  a  method  based  on  the  conser¬ 
vation  of  energy  and  momentum.  This  method  overcomes  the 
aforementioned  difficulty. 
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1. 


Introduction 


In  this  paper  we  develop  a  method  for  the  numerical  solution 
of  the  equations  of  motion  of  an  object  acted  upon  by  several 
gravitational  masses.  In  general,  the  motion  can  be  described  by 
a  special  class  of  second  order  differential  equations,  namely 

y"(x)  =  f(x,y(x))  .  (1) 

There  exist  methods  of  Runge-Kutta  type  which  tackle  this 
problem  directly  (Collatz  [8,  p.  61],  de  Vogelaere  [23],  Scraton 
[20]  and  others).  Also,  linear  multistep  methods  of  the  form 

k  k 1 

Z  ajYn+j  =  h2  "  Bjfn+j  (2) 
j=0  j=0 

exist.  See,  for  example,  Henrici  [11,  p.  289],  Lambert  [13,  p. 
252]  and  others.  One  of  the  authors  has  developed  hybrid  methods 
for  such  classes  (Neta  and  Lee  [15]  ,  Neta  [16,  18] ) . 

The  direct  application  of  methods  of  class  (2)  to  problem 
(1),  rather  than  the  application  of  a  conventional  linear 
multistep  method  to  an  equivalent  first-order  system  is  usually 
recommended  (Ash  [1]). 

Special  methods  based  on  a-priori  knowledge  of  the  frequency 
were  developed  by  Bettis  [3],  Steifel  and  Bettis  [22],  Gautschi 
[10],  Neta  and  Ford  [17],  Neta  [19],  van  der  Houwen  and  Sommeijer 
[12],  Lyche  [14]  and  Sommeijer  et  al  [21]. 
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To  avoid  the  use  of  the  frequency  other  methods  were 
developed  satisfying  the  so-called  P-stability.  This  idea  of  P- 
stability  was  introduced  by  Lambert  and  Watson  [13].  They  showed 
that  such  methods  are  necessarily  implicit  and  cannot  have  order 
greater  than  two.  Chawla  [5]  and  Cash  [4]  independently  showed 
that  this  order-barrier  can  be  crossed  over  by  considering  hybrid 
two-step  methods.  Other  P-stable  methods  were  developed  by 
Costabile  &  Costabile  [9],  Chawla  [6],  Chawla  and  Rao  [7]  and 
others . 

Our  idea  here  is  to  develop  a  new  method  to  approximate  the 
solution  of  the  two  body  problem.  This  scheme  will  conserve  both 
the  energy  per  unit  mass  and  the  specific  angular  momentum  of  the 
system. 

In  the  next  section  the  method  is  developed  for  the 
perturbation  free  flight.  Numerical  experiments  with  the  scheme 
will  be  described  in  the  last  section. 


2.  Development  of  the  Method 

In  this  section  we  describe  our  new  scheme  to  solve  the 
following  system  of  equations 


”  2  2, 3/2 

(x  +y  ) 


(3) 


,  2  2,3/2 

(x  +y  ) 


(4) 


where  k  =  3.9877  1014  m3  sec-2  is  the  gravitational  parameter. 

It  is  well  known  that  Cowell's  method  gives  a  numerical 
solution  that  spirals  inward.  Encke's  method  improves  this 
result  but  requires  more  work  [2]. 

It  can  be  easily  shown  (see  e.g.,  Bate  et  al.  [2])  that  the 
energy  and  momentum  are  conserved  for  a  perturbation  free  flight. 
The  conservation  of  these  quantities  will  be  used  to  obtain  the 
fast  changing  variable.  It  is  thus  important  to  rewrite  the 
system  in  polar  coordinates 


,d6. 2 

r(dt) 


k_ 

2 


(5) 


d^9  _  2_  dr  d£  (6) 

,  .2  r  dt  dt 

dt 


The  radius  r  does  not  change  much  in  time  and  thus  we  can  use  any 
method;  here  we  use  a  second  order  Taylor  series  method.  The 
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value  of  6  corresponding  to  r  will  be  obtained  by  integrating  the 
relation  (conservation  of  energy  per  unit  mass) 


(d£)2  +  (r  — ) 2 

ldt'  1  dt 


<afl  /  + 

t  =  0 


lr,0>  H't.0* 


To  be  more  specific,  our  method  approximates  r  using 


r  (ti+i) 


=  r(ti)  +  r(ti) (At) 


rr  (t .  )  (At) 


where 


r(ti+l)  =  rfti)  +  r(i)(At)  (9) 

and  r(t)  is  obtained  from  the  differential  equation.  In  a 
similar  fashion  we  find 

e(ti+l}  =  5(ti)  +  5 ( ti ) ( At )  +^9(ti)(At)2  (10) 

where 

S(t.)  =  \[7c  +  k/r(ti))2  -  r2(ti)/r(ti)  (11) 

and  6 ( t i )  is  obtained  by  differentiating  the  above  relation. 

The  results  of  this  method  were  compared  to  Taylor  series 
method  for  both  r  and  9  where  now  0  is  evaluated  from  an  equation 
similar  to  (9) . 
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any 


Bsmaris:  It  should  be  clear  that  this  idea  can  be  used  for 
conservative  system.  One  uses  the  conservation  of  some  quantity 
to  obtain  an  approximate  value  for  the  faster  changing  component. 
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3.  Numerical  Experiments 

In  this  section  we  solve  the  equations  of  motion  (3) -(4) 
combined  with  the  initial  values 

r ( 0)  =  R  +  105  , 

r(0)  =  0  ,  (12) 

6(0)  =  0  , 

5(0)  =  8000/r (0)  , 

where  R  =  6.371v106  is  the  radius  of  earth. 

In  the  first  table  we  compare  the  results  of  our  method  with 
Taylor  series  using  t  =  1  sec.  As  can  be  seen  in  the  table,  the 
value  of  r  at  the  end  of  each  period  is  nearly  constant  with  our 
method  but  decreases  as  the  number  of  periods  increases  with 
Taylor  series  method. 

TABLE  1 


At  the 

7 

Energy  per 

Specific  angul 

end  of 

r  no 

m) 

unit  mass 

momentum 

Period 

Taylor 

New 

(108  m^/sec^) 

(1011  m2/sec) 

0 

.  6471 

.  6471 

-.2962 

.  5277 

5 

.  6320 

.6476 

-.3019 

.5127 

10 

.6173 

.  6482 

-.3077 

.  5078 

15 

.6030 

.6487 

-.3135 

.  5029 

20 

.5890 

.  6493 

-.3194 

.4982 

30 

.  5622 

.6503 

-.3313 

.4888 

40 

.  5367 

.  6513 

-.3433 

.4798 

50 

.5123 

.  6523 

-.3557 

.4710 

60 

.4891 

.6533 

-.3683 

.  4624 

70 

.4667 

.6543 

-.3811 

.4540 
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The  energy  per  unit  mass  E  and  the  specific  angular  momentum  are 
given  by 


E 


1*2  *  2 

J(r^  +  (rep) 


k 

r  ' 


(13) 


P  = 


2  * 

re  -  r  9 


In  our  example  both  the  energy  and  momentum  are  conserved. 

in  the  following  three  figures  the  radius  r,  the  energy  per 
unit  mass  E  and  the  specific  angular  momentum  P  were  plotted  at 
the  end  of  every  10  periods.  It  is  clear  that  the  relative 
change  in  r  values  after  70  periods  using  Taylor  series  method  is 
more  than  25  times  larger  compared  to  our  new  method.  The  energy 
and  momentum  were  kept  constant  by  the  new  sc  .erne  whereas  the 
Taylor  series  method  shows  a  relative  loss  of  28%  in  energy  and 
12%  in  momentum  at  the  end  of  70  periods. 

Remark;  Note  that  our  method  is  explicit  whereas  the  P-stable 
methods  suggested  by  Lambert  and  Watson  [13]  are  implicit. 


COMPARISON  OF  RADIUS 


TIME  (PERIODS) 


COMPARISON  OF  MOMENTUM 


I  I MK  (PKKIODS) 
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